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ABSTRACT 

We use N-body/gasdynamical cosmological simulations to examine the effect of the 
assembly of a central galaxy on the shape and mass profile of its surrounding dark mat- 
ter halo. Two series of simulations are compared; one that follows only the evolution 
of the dark matter component of individual halos in the proper ACDM cosmological 
context, and a second series where a baryonic component is added and followed hy- 
drodynamically. The gasdynamical simulations include radiative cooling but neglect 
the formation of stars and their feedback. The efficient, unimpeded cooling that re- 
sults leads most baryons to collect at the halo center in a centrifugally-supported disk 
which, due to angular momentum losses, is too small and too massive when compared 
with typical spiral galaxies. This admittedly unrealistic model allows us, nevertheless, 
to gauge the maximum effect that galaxies may have in transforming their surrounding 
dark halos. We find, in agreement with earlier work, that the shape of the halo becomes 
more axisymmetric: post galaxy assembly, halos are transformed from triaxial into es- 
sentially oblate systems, with well-aligned isopotential contours of roughly constant 
flattening ((c/a) ~ 0.85). Halos always contract as a result of galaxy assembly, but 
the effect is substantially less pronounced than predicted by the traditional "adiabatic 
contraction" hypothesis. The reduced contraction helps to reconcile ACDM halos with 
constraints on the dark matter content inside the solar circle and should alleviate the 
long-standing difficulty of matching simultaneously the scaling properties of galaxy 
disks and the galaxy luminosity function. The halo contraction we report is also less 
pronounced than found in earlier simulations, a disagreement that suggests that halo 
contraction is not solely a function of the initial and final distribution of baryons. Not 
only how much baryonic mass has been deposited at the center of a halo matters, but 
also the mode of its deposition. It might prove impossible to predict the halo response 
to galaxy formation without a detailed understanding of a galaxy's detailed assembly 
history. 
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1 INTRODUCTION 

Over the past couple of decades, cosmological N-body sim- 
ulations have established a number of important results re- 
garding the structure of cold dark matter (CDM) halos. 
In particular, broad consensus holds regarding the overall 
shape of the halo mass distribution and its radial mass pro- 
file. CDM halos are distinctly triaxial systems supported by 
anisotropic velocity dispersion tensors (Frenk et al., 1988; 
Jing & Suto, 2002; Allgood et al., 2006; Hayashi et al., 2007; 
Bett et al., 2007), a result that reflects the highly aspher- 
ical nature of gravitationally-amplified density fluctuations 
as they become non linear and assemble into individual ha- 



los. Despite their chaotic assembly, the spherically-averaged 
mass profile of CDM halos is nearly "universal", and can 
be approximated fairly accurately, regardless of halo mass 
or redshift, by scaling a simple formula (see, e.g., Navarro 
et al., 1996b, 1997, 2004, 2008). 

One major limitation is that these results have been 
obtained from simulations that neglect the presence of the 
baryonic, luminous component of galaxy systems. Although 
this may be a suitable approximation in extremely dark 
matter-dominated systems, it is expected to fail in regions 
where baryons contribute a substantial fraction of the mass, 
as is the case in the luminous regions of normal galaxies 
like the Milky Way. The response of the dark halo to the 
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assembly of a galaxy is thus a crucial ingredient of models 
that attempt to interpret observational constraints in terms 
of the prevailing CDM paradigm. 

There has been extensive work, both numerical and an- 
alytical, on the modification of the mass profile of a dark 
halo induced by the assembly of a central galaxy. Early re- 
sults suggested that simple analytical prescriptions based 
on the conservation of adiabatic invariants gave an accu- 
rate description of the halo response. Following the early 
work by Barnes & White (1984), Blumenthal et al. (1986) 
devised a simple formula to link the dark mass profiles 
before and after the assembly of a galaxy. Given the ini- 
tial, spherically-symmetric enclosed mass profiles of the dark 
matter, M^ m (r), and baryons, M£(r), one may derive the fi- 
nal dark mass profile, M^ m (r), once the final baryonic mass 
profile, M[ (r), is specified. The model assumes that dark 
matter particles move on circular orbits before and after the 
contraction, and that their initial, r;, and final, 77, radii are 
related by the condition: 

H [Ml (n) + M dm ] = n [M dm + Mi{n)l (1) 

where Mdm = Afjmfr/) = Md m (7*i) is the dark mass en- 
closed by each dark matter particle (i.e., no shell crossing). 

Despite its simplicity, and the crude approximation on 
which it is based, early N-body work (e.g., Barnes & White, 
1984; Blumenthal et al., 1986; Jesseit et al., 2002) found 
reasonable agreement between eq. 1 and the results of simu- 
lations, and helped to establish the "adiabatic contraction" 
formulation as the default procedure when considering the 
halo response to the formation of a central galaxy. 

Semianalytic models of galaxy formation have adopted 
this formulation in order to link the (observed) dynamical 
properties of luminous galaxies to the (predicted) properties 
of the dark halos that surround them (see, e.g., Mo et al., 
1998; Cole et al., 2000; Dutton et al., 2005; Croton et al., 
2006; Somerville et al., 2008). This work has highlighted a 
number of potential problems when attempting to reconcile 
the predicted properties of galaxies in the ACDM cosmogony 
with observed scaling relations. 

One particularly important challenge has been the in- 
ability of semianalytic galaxy formation models to match 
simultaneously the zero point of the Tully-Fisher (TF) rela- 
tion and the galaxy luminosity function. Successful models 
require that the rotation speed of disks be of the order of the 
virial velocity of the halos they inhabit (see, e.g., Somerville 
& Primack, 1999; Croton et al., 2006). However, models that 
include adiabatic contraction typically predict disk rotation 
speeds well in excess of the virial velocity. 

As discussed most recently by, for example, Dutton 
et al. (2005); Gnedin et al. (2007) and Dutton et al. (2008), 
potential solutions to this problem include: (i) revising the 
adiabatic hypothesis so as to reduce (or even reverse!) the 
halo contraction; (ii) adopting lighter stellar mass-to-light 
ratios in order to minimize the contribution of baryons and 
to allow for more dark mass enclosed within disk galaxies; or 
(iii) modifying the cosmological parameters so as to reduce 
halo concentrations. We shall concentrate our analysis here 
on option (i), although we note that all three possibilities 
are probably equally important and should be treated on 
equal footing when addressing these issues in semianalytic 
models of galaxy formation. 

The possibility that the actual response of the halo to 



the assembly of the disk might differ substantially from the 
predictions of the adiabatic contraction formalism has been 
noted before. Barnes (1987) and Sellwood (1999) have re- 
marked that the adiabatic contraction hypothesis might lead 
to substantial overestimation of the halo compression. As 
discussed by Sellwood & McGaugh (2005) and Choi et al. 
(2006), such deviations are likely to depend strongly on the 
orbital structure of a halo, performing best when most par- 
ticles have large tangential motions, but poorly in systems 
with radially anisotropic velocity distributions. 

It should also be pointed out that the studies men- 
tioned above were carried out by perturbing simple spher- 
ical models with a potential term designed to imitate the 
growth of an assembling galaxy. These studies, therefore, 
miss the hierarchical nature of the assembly of a galaxy and 
of its surrounding halo. The study of Gnedin et al. (2004), 
on the other hand, uses the proper cosmological context, 
but focusses on simulations of galaxy clusters, rather than 
the galaxy-sized systems of interest for the issues discussed 
above. In this respect, our study is similar to that of Gustafs- 
son et al. (2006), who studied four simulations of galaxy- 
sized halos in the ACDM scenario. 

The assembly of a central galaxy also modifies the three- 
dimensional shape of the surrounding dark matter halo. This 
was already noted in the first simulations to include, in ad- 
dition to a dark matter halo, a dissipative baryonic compo- 
nent, such as the early work of Katz & Gunn (1991) and 
Katz & White (1993). Dubinski (1994) studied this further 
by growing adiabatically a central mass concentration inside 
a triaxial dark matter halo and confirmed that the steepen- 
ing of the potential leads to much rounder halo shapes. This 
result was confirmed by Kazantzidis et al (2004), who also 
noticed the effect in their cluster simulations. We extend this 
body of work by focussing, as in Hayashi et al. (2007), on 
the shape of the gravitational potential rather than on the 
axial ratios of the inertia tensor, as well as on its dependence 
on radius. 

The plan of this paper is as follows: § 2 describes the 
numerical simulations and § 3 presents the main results. 
After a brief general description of the evolution in § 3.1, 
§ 3.2 describes the main properties of the simulated central 
galaxies; § 3.3 discusses the modifications they induce on 
the shape of the halo gravitational potential; while § 3.4 
compares the mass profiles of dark halos before and after 
the inclusion of a dissipative baryonic component. We apply 
these results to the Milky Way in § 3.5 and compare with 
earlier work in § 3.6. We end with a brief summary of our 
main conclusions in § 4. 



2 NUMERICAL SIMULATIONS 
2.1 Cosmological parameters 

We adopt cosmological parameters consistent with the com- 
bined analysis of the 2dfGRS (Colless et al., 2001) and the 
first-year WMAP data (Spergel et al., 2003): a present day 
value of the Hubble constant of Ho = 70 km/s/Mpc; a scale- 
free initial density perturbation spectrum with no tilt and 
normalized by the linear rms mass fluctuations on 8fe _1 
Mpc spheres, as = 0.9. The matter-energy content of the 
Universe is expressed in units of the present-day critical 
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Figure 1. Dark matter (black) and gas (colored) particles for halo S02h. The color of a gas particle reflects its temperature, T < 10 4,5 K 
(green) or T > 10 45 K (magenta). Each panel zooms into the center of the system by consecutive factors of 3 in radius (see axis labels). 
The circle in the top-left panel shows the virial radius. The circle in the bottom-right panel shows the radius, r g i x , used to define the 
central galaxy. Note that the colder (green) gas component inhabits the center of the main halo and of its substructures, where it forms 
thin, centrifugally-supported disks. The hotter component (magenta) is distributed more or less uniformly across the main halo and 
makes up only 16% of all the baryons within the virial radius. 



density for closure, and contains a dominant cosmological 
constant term, Oa = 0.7, as well as contributions to the 
matter content, Qm = 0.3, from cold dark matter (CDM), 
ficDM = 0.259, and baryons, Sib = 0.041. 

2.2 Code and initial conditions 

We have performed a suite of 13 numerical simulations of the 
formation of galaxy-sized halos in the ACDM cosmogony. 
Each simulation follows the evolution of a relatively small 
region of the Universe, excised from a 432 3 -particle simula- 
tion of a large 50 /i -1 Mpc periodic box (Reed et al., 2003), 
and resimulated at much higher mass and spatial resolution. 



Each "zoomed-in" re-simulation follows the formation of a 
single galaxy-sized halo of mass ~ 10 12 Mq and its immedi- 
ate surroundings. The simulations include the tidal fields of 
the parent simulation, and follows the coupled evolution of 
gas and dark matter. The hydrodynamical evolution of the 
gaseous component is followed using the Smooth Particle 
Hydrodynamics (SPH) technique. This re-simulation tech- 
nique follows standard practice, as described, for example, 
by Power et al. (2003). All our simulations were performed 
using GASOLINE, a parallel N-body/SPH code described 
in detail by Wadsley et al. (2004). 



© 2008 RAS, MNRAS 000, 000-000 



4 Abadi et al 




Figure 2. Radius containing 90% of the baryonic mass, rgo, ver- 
sus the total baryonic mass of the central galaxy, M^k, in our 
simulated halos (filled circles), expressed in units of the halo virial 
values. The dashed vertical line the universal baryon fraction, 
/bar = ^b/^M- About 70-80% of all baryons have been collected 
in the central galaxy. The solid curve shows the corresponding val- 
ues for the Milky Way, varying the virial mass of its halo. Symbols 
along the line indicate a few different values of the virial veloc- 
ity, V vir = 220, 150, 110, and 91 km/s. Note that the simulated 
central galaxies are substantially more massive and smaller than 
spirals like the Milky Way. This maximizes the effect of baryons 
in transforming the shape and mass profile of the dark halo. 



Figure 3. Angular momentum and mass of the central galaxy in 
our simulations, expressed in units of the virial values. The curve 
labeled Jdisk = jvii corresponds to central galaxies with the same 
specific angular momentum as the surrounding halo (a common 
assumption of scmianalytic models). Those on the curve labeled 
jdisk/ivir = Tt<i//bar have specific angular momenta (in units of 
the halo's) that scale in proportion to the fraction of baryons 
within the virial radius that have collected in the central galaxy. 
Note that the simulated central galaxies have even lower angular 
momenta than in the latter assumption, highlighting the large 
angular momentum losses that accompany the assembly of the 
central galaxy, and explaining the small sizes of the disks shown 
in Fig. 2. 



2.3 Halo selection and resimulation 

The 13 halos were selected from a list of all ~ 10 12 Mq 
halos in the parent simulation, with a mild bias to avoid 
objects that have undergone a major merger after z — 1 or 
that have, at z = 0, unusually low (A < 0.03) spin param- 
eter. The mass resolution of the resimulations is such that 
halos are represented with ~ 40, 000 dark matter particles 
within their virial radius * at z = 0. Two sets of simu- 
lations are performed for each halo, one where only a dark 
matter component is followed, and another where dark mat- 
ter and baryons are included. Dark matter-only simulations 
("DMO" for short) assume that the total matter content of 
the Universe is in cold dark matter, Qm = ^cdm = 0.3. 

Simulations including a baryonic component ( "DM+B" 
for short) include equal number of gas and dark matter parti- 
cles, with masses modified to satisfy our choice of Q b = 0.041 
and ficDM = 0.259. The gaseous component is evolved us- 

* We define the virial radius, r v ; r , of a system as the radius of 
a sphere of mean density A v j r (,z) times the critical density for 
closure. This definition defines implicitly the virial mass, M v i r , as 
that enclosed within r v ; r , and the virial velocity, V v i T , as the circu- 
lar velocity measured at r v i r . Quantities characterizing a system 
will be measured within r v ; r , unless otherwise specified. The virial 
density contrast, A v j r (z) is given by A v i T (z) = 18n 2 + 82/(z) — 
39/(z) 2 , where f(z) = [Q (l + z) 3 /(n (l + z) 3 + ^a))] - 1 and 
Hq = !i C DM + Of, (Bryan & Norman, 1998). A vir 100 at z = 0. 



ing SPH including, besides the self-gravity of gas and dark 
matter, pressure gradients and shocks. Baryons are allowed 
to cool radiatively according to the cooling function of a gas 
of primordial composition down to a temperature of 10 4 K, 
below which cooling is disabled. No star formation or feed- 
back is included. As we discuss below, this choice maximizes 
cooling and favors the collection of most baryons at the cen- 
ter of each dark halo. Although unrealistic as a model for 
galaxy formation, this choice has the virtue of simplicity and 
also allows us to examine the halo response when the effect 
of the baryonic component is maximal. 

Pairwise gravitational interactions are softened adopt- 
ing a spline softening length kept fixed in comoving coordi- 
nates. We test for numerical resolution effects by simulating 
each halo with 8 times fewer particles. In the interest of 
brevity we do not present results from this low-resolution 
series here, but we have checked explicitly that none of the 
conclusions we present here are modified by this change in 
numerical resolution. We have also increased the number of 
particles in one case. Halo S02h is equivalent to S02, but 
with ~ 3.4 times more particles: at z = this halo has 
~ 140, 000 particles per component within the virial radius. 

Table 1 lists the main properties of the simulated halos. 
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Figure 4. Dark matter particles corresponding, at z = 0, to the two resimulations of halo S02h, "dark matter only" (left panels) 
and "dark matter plus baryons" (right panels). Particles are colored according to the value of the gravitational potential (binned in 
logarithmic units, top panels) or the local density (bottom panels). Note that the halo responds to the assembly of the central galaxy 
(shown in black in the right-hand panels) by becoming noticeably more spherical. Projections are chosen so that the rotation axis of the 
central disk coincides with the z-axis. The central disk is shown with black dots, to illustrate the orientation of the disk relative to the 
halo shape. In general, the central disk is well aligned with the minor axis of the halo. 



3 RESULTS 

3.1 General evolution 

Because our simulations neglect the formation of stars and 
their feedback, the evolution of the baryonic component is 
characterized by the rapid cooling and collapse of baryons 



at the center of the early collapsing progenitors of the fi- 
nal halo. As a result, the main mode of galaxy assembly is 
mergers: up to 70% of all baryons in the central galaxy were 
accreted in the form of dense, cold, gaseous clumps that sink 
to the center through dynamical friction. 

As shown in earlier work, (see, e.g., Navarro & Benz, 
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Figure 5. "Radial" (r = V a' 2 + b 2 + c 2 ) dependence of axial 
ratios b/a (open symbols) and c/a (filled symbols) for our high- 
resolution run S02h. Halo shapes are measured by fitting 3D el- 
lipsoids to the position of dark matter particles in narrow loga- 
rithmic bins of the gravitational potential. The potential is com- 
puted using only dark matter particles, so it actually corresponds 
to the contribution of the dark matter component to the overall 
potential. Red (lower) symbols correspond to the DMO run; blue 
(upper) symbols to the DM+B run. Note that the central galaxy 
turns a rather triaxial, nearly prolate halo into an axisymmetric, 
nearly oblate one. Curves are fits using eq. 2. Fit parameters for 
this halo and the mean value computed over the sample are listed 
in Table 2. 




0.5 1 1.5 2 2.50 0.5 1 1.5 2 2.5 
Log(r/kpc) Log(r/kpc) 



Figure 6. "Radial" (r = \/d 2 + b' 2 + c 2 ) dependence of axial 
ratios b/a (upper panel) and c/a (lower panel) for all simulated 
haloes. The curves for the "dark matter only" runs are computed 
by fitting eq. 2 to the radial dependence of the axial ratios of par- 
ticles along isopotential contours. The contours in all cases refer 
to the potential contributed solely by the dark matter compo- 
nent. Note that, as a result of the assembly of the central galaxy, 
(i) halos become nearly oblate and (ii) axial ratios become ap- 
proximately independent of radius. The effect on the shape of the 
halo extends almost to the virial radius, far beyond the actual 
size of the central galaxy. Thick lines correspond to mean values 
computed over the sample. 



1991; Navarro & White, 1994; Navarro et al., 1997), this 
mode of assembly leads to very efficient accretion of baryons 
into the central galaxy and to the transfer of much of 
their angular momentum to the surrounding halo. Because 
baryons are assumed to remain in gaseous form at all times, 
those that can cool are forced to settle into centrifugally- 
supported disks at the center of their surrounding dark ha- 
los. Thus, we expect the central galaxies in our simulations 
to be disks that contain a large fraction of all baryons within 
the virial radius, and to have angular momenta well below 
the angular momentum content of the system as a whole. 

Figure 1 illustrates this for the case of our highest reso- 
lution halo, S02h (see Table 1). The four panels in this figure 
show the dark matter (black) and baryonic (colored) compo- 
nents, zooming by consecutive factors of 3 toward the central 
galaxy. Cold gas (T < 10 4,5 K) is shown in green, whereas 
hot gas (T > 10 4 ' 5 K) is shown in magenta. Note that most 
of the cold gas inhabits the center of the main halo and its 
substructures, where it forms easily identifiable thin, mas- 
sive disks. All cold baryons associated with the central disk 
are contained within a sphere of radius r g i x = 10 kpc, which 
we shall use hereafter to define the central galaxy in all runs. 



3.2 Mass, size and angular momentum of central 
galaxies 

The central galaxy contains most of the baryons within the 
virial radius of the system (74% in the case of S02h, see 
Table 1). This is shown in Fig. 2, where we plot, for all our 
simulations, the baryonic mass of the central galaxy, A/disk, 
versus the radius, rgo, that contains 90% of its mass, scaled 
to the virial mass and radius of the system, respectively. The 
vertical dashed line indicates the universal baryon fraction 
in the simulations, /bar = fibar/fiM = 0.137. Typically, 70 
to 80% of all baryons within r v i r are found in fairly small 
central disks that are fully contained within a radius of order 
~ 3% of the virial radius. 

We may compare this with a typical spiral galaxy like 
the Milky Way (MW), where the mass and size of the bary- 
onic component can be estimated accurately. Assuming that 
the disk is exponential with total mass 4.5 x 10 10 Mq and 
radial scalelength R^f w — 2.5 kpc and that the mass of the 
bulge is A'/buigc = 4.5 x 10 9 Mq, we estimate that 90% of the 
Milky Way baryons are confined within rgo W = 9.4 kpc. 

In order to compare this with the simulation results 
shown in Fig. 2 we need make an assumption about the 
virial mass of the Milky Way halo. For a virial velocity of 
the order of the rotation speed at the solar circle (220 km/s), 
the baryonic component makes up only ~ 1% of the virial 
mass, or roughly 7% of all baryons within r v i r . The simulated 
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disks are thus much smaller than a typical spiral like the 
MW. They are also comparatively more massive. 

Only if the virial velocity is as low as 110 km/s do the 
MW bulge+disk make up a fraction of available baryons 
as high as in our simulations (~ 75%). However, in that 
case r^J w ~ 0.06r v i r , a factor of ~ 4 times larger than our 
simulated disks. 

The reason why simulated gaseous disks are so much 
smaller than typical spirals is that baryons have lost a large 
fraction of their angular momentum to the surrounding halo. 
This is illustrated in Fig 3. where we plot the angular mo- 
mentum of the gaseous disk (in units of that of the system 
as a whole, which is dominated by the dark matter com- 
ponent), J d = Jdisk/Jvir = M disk j disk /M vil j vil , versus the 
mass fraction, rtid = Mdi S k/M v i r . These are the parameters 
commonly adopted in semianalytic models of disk formation, 
such as those of Mo et al. (1998). 

Since baryons acquire during the expansion phase as 
much angular momentum as the dark matter, and it is un- 
likely that Mdisk/Mvir will exceed the universal baryon frac- 
tion, then we do not expect the specific angular momentum 
of the disk, jdisk, to exceed that of the system as a whole, 
jvir- Therefore, central galaxies are unlikely to populate the 
shaded areas of Fig. 3. 

Most semianalytic work assumes, for simplicity, that 
jdisk = jvir or, equivalently, that J d — m d , regardless of the 
value of md- On the other hand, Navarro & Steinmetz (2000) 
argue that it is unlikely that a galaxy where rrid < < /bar may 
have the same specific angular momentum as the whole sys- 
tem, and propose that jdisk, as a fraction of j v j r , should be 
comparable to the fraction of baryons that make up the cen- 
tral galaxy; i.e., jdisk/jvir = M diBk /(/ bar M vir ) = m d // bar . 
This is illustrated in Fig. 3 by the lower diagonal line. 

Our simulated galaxies are well below both lines. This 
indicates that baryons have specific angular momenta much 
lower than their surrounding halos and explains their small 
sizes compared to typical spirals. These are therefore unreal- 
istic models of spiral galaxy formation, but the combination 
of large mass and small size allows us to probe the response 
of the dark halo in the case where the deepening of the po- 
tential well due to the central galaxy is maximal. 

3.3 Halo shape 

As anticipated in § 1, the halo responds to the presence of 
the central galaxy by becoming significantly more spherical. 
This is illustrated in Fig. 4, where we compare the shape 
of isodensity and isopotential contours for the "dark matter 
only" and "dark matter plus baryons" runs of halo S02h. 
Panels on the left correspond to the DMO run, those on 
the right to the DM+B run. Particles are colored according 
to their local values of their gravitational potential or local 
density (binned in a logarithmic scale), computed using only 
the dark matter particles. Note that using the gravitational 
potential leads to much more stable estimates of the shape 
of the halo, since it is much less affected by the presence of 
substructures and other transient fluctuations in the mass 
distribution. 

The isopotential contours are well approximated by el- 
lipsoidal surfaces, and we use the axial ratios of such el- 
lipsoids to measure, as a function or radius, the change in 
shape of halo S02h. We show the result in Fig. 5. As dis- 



cussed by Hayashi et al. (2007), the radial dependence of 
the axial ratios may be approximated by the formula, 



log(- or -) 

a a 



tanh 



( 7l0g ^) 



1 



(2) 



Here a parameterizes the central value of the axial ratio, 
(b/a)o or (c/a)o, by 10~ 2 a ; r a indicates the characteristic 
radius at which the axial ratio increases significantly from its 
central value; and 7 regulates the sharpness of the transition. 

The presence of the central galaxy turns the halo from a 
triaxial, nearly prolate system into an axisymmetric, nearly 
oblate system where the axial ratio is nearly constant with 
radius. As a result, eq. 2 is not adequate for the nearly fea- 
tureless radial dependence of the axial ratios in the DM+B 
runs. We fit the latter with a simple power-law, 

b c 

log(— or — ) = alogr + j3, (3) 
a a 

The best-fit parameter to the mass and isopotential con- 
tour profile shapes are given in Table 2. 

The profile fits for all simulations are compiled in Fig. 6, 
and illustrate a few important points. As a result of the as- 
sembly of the central galaxy: (i) halos become nearly oblate; 
(ii) axial ratios are roughly independent of radius; and (iii) 
halo shapes are affected well beyond the size of the central 
galaxy, and nearly as far out as the virial radius. 



3.4 Mass profile and contraction 

Fig. 7 shows the enclosed mass profile of the DMO and 
DM+B runs corresponding to halo S02h. DMO dark masses 
have been scaled by (1 — /bar) so that the total dark mass 
in both the DMO and DM+B runs are the same. The DMO 
mass profile is shown by the thick red solid curve and, as ex- 
pected, it is well approximated by the NFW (Navarro et al., 
1996b, 1997) formula (thin solid line). The central disk that 
assembles in the DM+B run leads to a contraction of the 
dark halo: there is more dark mass in the DM+B run at all 
radii compared with the DMO run. This is a result com- 
mon to all our simulations; in no case do we see the halo 
"expand" as a results of the assembly of the central galaxy. 

The contraction, however, is not as pronounced as what 
would be expected from the "adiabatic contraction" model 
discussed in § 1. The adiabatic contraction prediction (eq. 1) 
is shown by the thick dotted line in Fig. 7. The discrepancy 
between model and numerical results is not small. At ~ 5 
kpc, a radius roughly twice the gravitational softening and 
that encloses more than 1300 dark particles, the adiabatic- 
contraction formula predicts ~ 2.5 times more dark mass 
than found in our simulations. Even at a radius of 10 kpc, 
eq. 1 overestimates the halo contraction by more than ~ 50% 
in mass. The modified contraction proposed by Gnedin et al. 
(2004) (see curve labeled "Contra" in Fig. 7) fares better, 
but it still overpredicts the results of the simulations. 

A similar result applies to all of our simulated halos. 
We illustrate this in Fig. 8, where we plot the ratio between 
radii that contain a given number of dark matter particles 
in the DMO run (n) and the DM+B run (77) versus the 
ratio between M; , the total mass within n in the DMO run, 
and Mf, the total mass within r; in the DM+B run. With 
this choice, the adiabatic-contraction prediction is simply 
r f/ r i — Mi/Mf, which is traced by the 1:1 line in Fig. 8. 



© 2008 RAS, MNRAS 000, 000-000 



8 Abadi et al 



12 -// 




0.5 1 1.5 2 2.5 

Log[r/kpc] 



Figure 7. Enclosed mass profile of various components of halo 
S02h. The "dark matter only" profile is shown with a thick curve, 
after scaling masses by (1 — /bar): so that, within r v ; r , the total 
dark mass of the DMO and DM+B runs will be comparable. The 
thin line shows a Navarro-Frenk- White halo fit to the DMO pro- 
file. Other colors and line types correspond to the various compo- 
nents of the DM+B run, as specified in the figure labels. The thick 
dashed blue curve shows the dark mass profile for the DM+B run. 
The assembly of most baryons into a central galaxy (dotted ma- 
genta curve) has clearly led to a contraction of the dark mass 
profile. The profile predicted by the "adiabatic contraction" for- 
mula (eq. 1; dot-dashed green curve) overpredicts the response 
of the halo. The modified adiabatic contraction model of Gnedin 
et al. (2004), shown by the dotted cyan curve, also overestimates 
the halo response. 



The numerical simulations are shown by the upper 
curves in this figure. Near the center, the baryons dominate 
over the DMO mass profile, and therefore Mi/Mf <C 1. Fur- 
ther out, the contribution of baryons to the total enclosed 
mass decreases, approaching the universal baryon fraction at 
the virial radius, where Mi/Mf tends to unity. The inner- 
most radius plotted in each case corresponds to that enclos- 
ing 1000 dark matter particles. This is in all cases comfort- 
ably larger than the gravitational softening, minimizing the 
possibility that our results are unduly influenced by numer- 
ical artifact. The open circles with error bars in Fig. 8 trace 
the median and rms scatter of all our simulation results. 

The ratio r//rj tends to a constant for Mi/Mf <C 1, 
suggesting that the halo response approaches saturation in 
regions where baryons dominate. A simple formula captures 
the average behaviour well (thick upper solid line), 

r / /r, = l + a[(M,/M f )"-l] 1 (4) 

with a = 0.3 and n = 2. As an application, we shall use this 
expression below to explore what constraints this implies for 
the mass and concentration of the halo of the Milky Way. 




0.2 0.4 0.6 0.8 



(M^r,)+M DM (r i ))/(M^r f )+M DM (r [ )) 

Figure 8. Dark halo response to the assembly of a central galaxy. 
The ordinate shows the ratio, rr/ri, between the radius contain- 
ing a given amount of dark mass in the DMO and DM+B runs, 
respectively, after scaling DMO masses by 1 — /bar- The smaller 
Tf/ri the stronger the halo contraction. The x-axis shows the ratio 
between the total mass, M;, contained within r, (in the DMO run) 
and Mf, that enclosed within r* (in the DM+B run). The adia- 
batic contraction formula (eq. 1) predicts that rf/r% = Mi/Mf, 
this is shown by the 1:1 line in the figure. Numerical results are 
shown for individual halos; from the radius that contains 1000 
dark particles outwards in order to minimize numerical uncer- 
tainties. The adiabatic-contraction formula overestimates the halo 
response. Gnedin et al. (2004) 's modified adiabatic contraction 
("Contra") does better but still overpredicts the halo response 
at most radii, and especially near the center. Symbols with error 
bars trace the median and quartiles of the numerical results. The 
thick upper curve is a fit using eq. 4. 



3.5 Application to the Milky Way 

The Milky Way offers a case study for the results described 
above. Because the mass and radial distribution of the bary- 
onic component, as well as the rotation speed of the local 
standard of rest (LSR), are relatively well known, the total 
dark mass contained within the solar circle is firmly con- 
strained. Adopting the same quantities for the bulge and 
disk of the Milky Way adopted in § 3.2, we find that baryons 
contribute (in quadrature) ~ 171 km/s to the circular ve- 
locity of the LSR, which we assume to be 220 km/s. This 
implies that the dark mass of the Milky Way within the solar 
circle (Rq = 8 kpc) is ~ 3.5 x 10 10 Mq. In order to allow for 
the possibility that some of this dark mass may be baryonic, 
we shall treat this a formal upper limit in the analysis that 
follows. 

The second constraint comes from the total baryonic 
mass of the Milky Way, under the plausible assumption that 
the mass of the central galaxy cannot exceed the total mass 
of baryons within the virial radius, ~ /barAf vlr . Thus, the 
minimum virial mass allowed for the MW halo by this con- 
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straint is 3.64 x 10 11 Mq, which corresponds to a virial ve- 
locity of ~ 92 km/s. 

Are these constraints compatible with ACDM halos? 
The top left panel of Fig. 9 shows the dark mass ex- 
pected within 8 kpc vs halo virial velocity. Each dot in this 
panel correspond to an NFW halo with concentration drawn 
at random from the mass-concentration relation (including 
scatter) derived by Neto et al. (2007) from the Millennium 
Simulation (Springel et al., 2005). As expected, the mass 
within 8 kpc increases with the virial velocity (mass) of the 
halo, modulated by the fairly large scatter in concentration 
at given halo mass. 

The black "wedge" in this panel indicates the region al- 
lowed by the MW constraints discussed above. As discussed 
by Eke et al. (2001), most ACDM halos satisfy the con- 
straints, and would be consistent with the MW if they were 
somehow able to avoid contraction. Note as well that the 
larger the MW halo virial velocity the lower the allowed con- 
centration. For V v ir — 220 km/s (the value required by semi- 
analytic models to match the Tully-Fisher relation and the 
luminosity function, see § 1) only halos with c v i r < 9.5 would 
be consistent with the Milky Way. For comparison, the av- 
erage concentration of Vvir = 220 km/s halos is (c v ir) ~ 10.4 
and its (lognormal) dispersion is ai ogc = 2.9, so this con- 
dition effectively excludes only 57% of such halos from the 
allowed pool. 

The situation is rather different if halos are adiabati- 
cally contracted (panel labeled "AdiabCont" in Fig. 9). The 
contraction increases substantially the dark mass contained 
within 8 kpc, so that very few halos satisfy the MW con- 
straints. For example, essentially no halo with Vvir ~ 220 
km/s could host a galaxy like the Milky Way, and the few 
that could would need to have c v ir < 3.2, many sigma away 
from the average concentration of halos of that mass. The 
halo of the Milky Way would need to be a very special halo of 
unusually low concentration if adiabatic contraction holds. 

Adopting the halo contraction of our numerical simula- 
tions improves matters. This may be seen in the bottom-left 
panel of Fig. 9, where we have contracted each halo using 
eq. 4. The range of allowed halo masses and concentrations 
is broader; even some halos with Vvir ~ 220 km/s could 
be consistent with the Milky Way, provided that c v i r < 6.5. 
This is significantly lower than the average at that mass, but 
conditions are much less restrictive if the halo of the Milky 
Way is less massive; for example, half of all Vvir ~ 130 km/s 
halos satisfy the MW constraints. The possibility that the 
virial velocity of the Milky Way is significantly lower than 
220 km/s has indeed been advocated by a number of recent 
observational studies (Smith et al., 2007; Sales et al., 2007; 
Xue et al., 2008). 

Lowering the average concentration at given mass also 
helps. This is shown in the bottom-right panel of Fig. 9, 
where we have repeated the exercise, but lowering all con- 
centrations by 20%. According to Duffy et al. (2008), this 
is approximately the change in average concentration when 
modifying the cosmological parameters from the values we 
adopt here (which are consistent with the first-year analy- 
sis of the WMAP satellite data) to those favoured by the 
latest 5-year WMAP data analysis. With this revision, 18% 
of Vvir ~ 220 km/s halos and roughly half of all Vvir ~ 160 
km/s halos would be consistent with the MW. 
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Figure 9. Dark mass within 8 kpc vs virial velocity for halos 
drawn at random from the ACDM mass-concentration relation 
(including scatter) of Neto et al. (2007). The DMO panel (top 
left) shows the results neglecting the effects of baryons and as- 
suming an NFW profile for the dark halos. The black "wedge" 
highlights the region of parameter space compatible with obser- 
vations of the Milky Way. The panel labeled "AdiabCont" shows 
the result of contracting each halo adopting the adiabatic contrac- 
tion formula (eq. 1). Very few halos would be consistent with the 
Milky Way if this formula holds. The bottom-left panels shows the 
result of applying our simulation results for the halo contraction, 
as captured by eq. 4. The bottom-right panel repeats the same 
exercise, but after reducing all concentrations by 20% in order to 
mimic the expected concentration change resulting from adopt- 
ing the latest cosmological parameters from the 5-year analysis 
of WMAP satellite data (Duffy et al., 2008). 



3.6 Comparison with earlier work 

The failure of the adiabatic contraction formalism to match 
the results of our simulations suggest that the response of 
the halo is more complex than what can be captured with a 
simple model for the contraction of spherical shells. This is 
confirmed by comparing our results with earlier work. 

Our results disagree not only with the traditional adi- 
abatic contraction formula, but also with the modified for- 
malism proposed by Gnedin et al. (2004) (shown with thin 
dash-dotted green lines in Fig. 8). These authors calibrated 
their results with numerical simulations not dissimilar to 
ours, except for the mass scale — they used mainly galaxy 
cluster halos. The mass scale, on the other hand, does not 
seem to be the reason for the discrepancy. Our halo contrac- 
tion is also less pronounced than found in two galaxy-sized 
halo simulations, as shown by the lines labeled A03&M03 
in Fig. 8. These correspond to gasdynamical simulations of 
galaxy formation presented by Abadi et al. (2003a,b) and 
Meza et al. (2003), and differ from ours mainly in their in- 
clusion of star formation and feedback effects. Those two 
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simulations seem to agree better with the "Contra" predic- 
tions than with our simulation results. 

Barring numerical artifact, these results seem to sug- 
gest that the halo response does not depend solely on the 
initial and final distribution of baryons. One possible expla- 
nation is that not only how much baryonic mass has been 
deposited at the center of a halo matters, but also the mode 
of its deposition. It is certainly plausible that central galax- 
ies assembled through merging of dense subclumps may lead 
to different halo response than galaxies assembled through 
a smooth flow of baryons to the center. Mergers of bary- 
onic subsystems may in principle pump energy into the dark 
halo (through dynamical friction) , altering its central struc- 
ture and softening its contraction. This possibility has been 
argued before (see, e.g., El-Zant et al., 2001; Ma & Boylan- 
Kolchin, 2004; Mo & Mao, 2004), and seems to find favor 
in our simulations, where mergers between massive baryonic 
clumps are frequent and plentiful and halo contraction is less 
strong than reported in earlier work. 

Recent work has also speculated that mergers may lead 
to halo expansion and that this would help to reconcile the 
properties of disk galaxies with ACDM halos (Dutton et al., 
2007). Our halos always contract, and it seems safe to con- 
clude that our results reflect the maximum effect of mergers 
on the halo response. We conclude therefore that mergers 
alone are unlikely to result in halo expansion. If such expan- 
sion is truly needed to reconcile disk properties with ACDM 
halos, it should come as a consequence of other processes not 
considered here, such as feedback-driven winds that may re- 
move substantial fraction of baryons from the central galaxy 
(see, e.g., Navarro et al., 1996a; Babul & Ferguson, 1996). 

4 SUMMARY 

We have used a suite of cosmological N-body/gasdynamical 
simulations to examine the modifications to the dark halo 
structure that result from the assembly of a central galaxy. 
The formation of 13 ACDM halos is simulated twice, with 
and without a baryonic (gaseous) component. For simplic- 
ity, the gasdynamic simulations include radiative cooling but 
neglect star formation and feedback. This favors the forma- 
tion of massive central baryonic disks at the center of the 
early collapsing progenitors of the final halo. As these sys- 
tems merge, the gaseous clumps merge and re-form a disk 
at the center of the remnant. 

The merger process leads to the transfer of a large frac- 
tion of the angular momentum from the baryons to the halo. 
At z = 0, the simulated central galaxies are too massive and 
too small to be consistent with observed spiral galaxies. Al- 
though unrealistic as a disk galaxy formation model, these 
simulations allow us to probe the dark halo response in the 
interesting case where the deepening of the potential well 
resulting from the formation of the central galaxy is maxi- 
mized. Our main conclusions may be summarized as follows. 

• Dark halos become significantly more spherical as a 
result of the assembly of the central galaxy. The triaxial, 
nearly prolate systems that form in the absence of a baryonic 
component are transformed into essentially oblate systems 
with a roughly constant axial isopotential ratio (c/a) ~ 0.85. 

• Halos always contract in response to the formation of 
the central galaxy. The "adiabatic contraction" formalism 



overestimates the halo contraction in our simulations. The 
discrepancy increases toward the centre, where the effect of 
baryons is larger. A simple empirical formula (eq. 4) de- 
scribes our numerical results. 

• The halo contraction in our simulations is also less pro- 
nounced than found in earlier numerical work (see, e.g., 
Gnedin et al., 2004), and suggest that the response of a halo 
does not depend on the final mass and radial distribution of 
baryons in the central galaxy, but also on the mode of their 
assembly. 

• We apply these results to the Milky Way, where accu- 
rate estimates of the mass of baryons and dark matter inside 
the solar circle exist. These allow us to probe the range of 
halo virial mass and concentration consistent with the con- 
straints. Only halos of unusually low mass and concentra- 
tion would match such constraints if "adiabatic contraction" 
holds. 

• This restriction is less severe if one uses the halo con- 
traction reported here (eq. 4): although few ACDM halos of 
virial velocity ~ 220 km/s would be eligible hosts for the 
Milky Way, the situation improves for lower virial veloci- 
ties. Halos of average concentration with virial velocity as 
large as V vlI = 135 km/s, would be consistent with the MW 
constraints. 

The dark halo response to the formation of a galaxy 
seems inextricably linked to the full coupled evolution of 
baryons and dark matter and may vary from system to sys- 
tem. Progress in our understanding of the distribution of 
dark matter in a baryon-dominated system will thus likely 
develop in step with our understanding of the particular as- 
sembly history of each individual system. 
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Label 


M vil 




-^DMO 


N DM 




Wgas 


JDMO 


JDM 


Jgas 




[10 12 M Q ] 


[kpc] 


(< r vir ) 


(< r vir ) 


(< rVi r ) 


(< r glx) 


[kpc km s _1 ] 


[kpc km s _1 ] 


[kpc km s _1 ] 


SOI 


1.22 


277.82 


53009 


52684 


54732 


38260 


1847.7 


2162.4 


462.0 


S02 


0.94 


254.18 


40847 


40687 


39775 


30613 


1063.3 


1088.3 


332.8 


S02h 


0.95 


255.89 


139748 


139729 


140618 


104309 


1095.7 


1200.3 


359.1 


S03 


1.51 


298.24 


65729 


65376 


66416 


37445 


2593.9 


2883.4 


438.1 


S04 


0.85 


246.02 


37419 


36432 


39007 


29534 


1443.1 


1626.7 


508.8 


S05 


0.99 


258.64 


42564 


42334 


45270 


30266 


2505.0 


2543.2 


173.3 


S06 


1.22 


278.05 


53509 


52652 


55674 


30411 


2736.0 


3056.0 


289.4 


S07 


0.86 


246.95 


37350 


37305 


36524 


25743 


1332.2 


1385.1 


310.0 


S08 


0.84 


245.39 


35902 


36224 


38221 


26031 


1221.3 


1292.5 


277.6 


S09 


0.84 


245.23 


36619 


36259 


37441 


28264 


850.5 


929.5 


531.8 


S10 


0.97 


257.68 


42004 


42074 


43469 


34475 


1976.6 


2231.8 


430.7 


Sll 


0.88 


249.53 


40181 


38190 


39551 


23579 


1430.3 


1676.8 


289.6 


S12 


0.90 


251.25 


39715 


39341 


38140 


28738 


715.7 


764.0 


381.5 


S13 


1.08 


266.72 


46863 


46337 


50223 


36708 


1470.7 


1607.5 


261.3 



Table 1. Main parameters of simulated halos. M v ; r is the total (dark matter plus gas) mass inside the virial radius, r v - lr . A^dmo an d 
)Vom is the number of dark matter particles inside the virial radius r v ; r in the "dark matter only" and "dark matter plus baryons" 
simulations, respectively. Columns 6 and 7 list the number of gas particles, -/V gas , within the virial radius and within the radius used to 
define the central galaxy, r g i x = 10 kpc. The specific angular momentum of the dark matter within r v i r and of the baryons within r g i x 
is also listed. 





b/a (Dark Matter Only) 


c/a (Dak Matter Only) 


b/a (Dark Matter+Gas) 


c/a (Dark Matter+Gas) 




a 7 r a 


a 7 r a 


a (3 


a (3 


S02h 


0.113 1.767 58.31 


0.137 2.102 83.13 


-0.00429 -0.00508 


0.00185 -0.07127 


Sample Mean Values 


0.081 1.383 35.81 


0.103 1.427 68.06 


-0.00934 -0.01009 


0.00081 -0.07220 



Table 2. Best fitting parameters for isopotcntial contour profile shapes. Columns 1-3 (4-6) correspond to axial ratio b/a (c/a) for Dark 
Matter Only runs (eq. 2). Columns 7-8 (9-10) correspond to axial ratio b/a (c/a) for Dark Matter+Gas runs (eq. 3). We quote values 
for the high resolution run S02h shown in figure 5 and also the mean value for the sample shown as thick lines in figure 6. 



Croton D. J., Springcl V., White S. D. M., et al., 2006, MNRAS, 
365, 11 

Dubinski J., 1994, ApJ, 431, 617 

Duffy A. R., Schaye J., Kay S. T., Dalla Vecchia C, 2008, MN- 
RAS, 390, L64 

Dutton A. A., Courteau S., de Jong R., Carignan C, 2005, ApJ, 
619, 218 

Dutton A. A., van den Bosch F. C, Courteau S., 2008, ArXiv 
e-prints 

Dutton A. A., van den Bosch F. C, Dekel A., Courteau S., 2007, 
ApJ, 654, 27 

Eke V. R., Navarro J. F., Steinmetz M., 2001, ApJ, 554, 114 
El-Zant A., Shlosman I., Hoffman Y., 2001, ApJ, 560, 636 
Frenk C. S., White S. D. M., Davis M., Efstathiou G., 1988, ApJ, 
327, 507 

Gnedin O. Y., Kravtsov A. V., Klypin A. A., Nagai D., 2004, 
ApJ, 616, 16 

Gnedin O. Y., Weinberg D. H., Pizagno J., Prada F., Rix H.-W., 

2007, ApJ, 671, 1115 
Gustafsson M., Fairbairn M., Sommer-Larsen J., 2006, Phys Rev 

D, 74, 12, 123522 
Hayashi E., Navarro J. F., Springel V., 2007, MNRAS, 377, 50 
Jesseit R., Naab T., Burkert A., 2002, ApJL, 571, L89 
Jing Y. P., Suto Y., 2002, ApJ, 574, 538 
Katz N., Gunn J. E., 1991, ApJ, 377, 365 
Katz N., White S. D. M., 1993, ApJ, 412, 455 
Ma C.-P., Boylan-Kolchin M., 2004, Physical Review Letters, 93, 

2, 021301 

Meza A., Navarro J. F., Steinmetz M., Eke V. R., 2003, ApJ, 590, 
619 

Mo H. J., Mao S., 2004, MNRAS, 353, 829 



Mo H. J., Mao S., White S. D. M., 1998, MNRAS, 295, 319 
Navarro J. F., Benz W., 1991, ApJ, 380, 320 

Navarro J. F., Eke V. R., Frenk C. S., 1996a, MNRAS, 283, L72 
Navarro J. F., Frenk C. S., White S. D. M., 1996b, ApJ, 462, 563 
Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493 
Navarro J. F., Hayashi E., Power C, ct al., 2004, MNRAS, 349, 
1039 

Navarro J. F., Ludlow A., Springel V., et al., 2008, ArXiv e-prints 
Navarro J. F., Steinmetz M., 2000, ApJ, 538, 477 
Navarro J. F., White S. D. M., 1994, MNRAS, 267, 401 
Ncto A. F., Gao L., Bctt P., ct al., 2007, MNRAS, 381, 1450 
Power C, Navarro J. F., Jenkins A., et al., 2003, MNRAS, 338, 
14 

Reed D., Gardner J., Quinn T., ct al., 2003, MNRAS, 346, 565 
Sales L. V., Navarro J. F., Abadi M. G., Steinmetz M., 2007, 

MNRAS, 379, 1464 
Sellwood J. A., 1999, in Galaxy Dynamics - A Rutgers Sympo- 
sium, edited by D. R. Merritt, M. Valluri, J. A. Sellwood, vol. 
182 of Astronomical Society of the Pacific Conference Series, 
351-+ 

Sellwood J. A., McGaugh S. S., 2005, ApJ, 634, 70 
Smith M. C, Ruchti G. R., Helmi A., ct al., 2007, MNRAS, 379, 
755 

Somcrvillc R. S., Barden M., Rix H.-W., et al., 2008, ApJ, 672, 
776 

Somcrvillc R. S., Primack J. R., 1999, MNRAS, 310, 1087 
Spergel D. N., Verde L., Peiris H. V., et al., 2003, ApJS, 148, 175 
Springel V., White S. D. M., Jenkins A., et al., 2005, Nature, 435, 
629 

Wadsley J. W., Stadel J., Quinn T., 2004, New Astronomy, 9, 137 
Xuc X. X., Rix H. W., Zhao G., et al., 2008, ApJ, 684, 1143 



© 2008 RAS, MNRAS 000, 000-000 



